
# coding: utf-8

# ## <b>1、首先明确学习 机器学习 的动机</b>
# -------------------------------------------------
# 先放一张大神的照片<br>
# 
# ![詹姆斯·西蒙斯](https://img1.doubanio.com/view/photo/photo/public/p2365189738.jpg)<br>
# 
# 不知道他的同学，点这里：[詹姆斯·西蒙斯][1]
# [1]:http://baike.baidu.com/link?url=f7gg3CfSgV0FnVShzVofMK65MC3UU6kfUbK0_RupJ0ZHr4rG2UCfYWz5PEHjxqodrhf0O22RKKXJhi3NiCVZxa
# 
# 詹姆斯·西蒙斯在TED的对话中有提到：<br>
# Q：<b>机器学习</b>在这里扮演了怎样一个角色？<br>
# A：某种意义上，<b>我们做的就是机器学习</b>。你观察一大堆数据，模拟不同的预测方案，<br>
# 直到你越来越擅长于此。我们所做之事，不见得一定有自我反馈，<b>但确实有效</b>。
# <p>
# 视频地址：[与横扫华尔街数学家的珍贵对话][2]
# [2]:http://v.qq.com/boke/page/s/0/w/s0186b0yrpw.html
# <p>
# 　　然而现在有许多已实现的机器学习开源包可供我们调用，如sklearn，更高级的技术还有Hadoop、Spark等等。<br>
# 那么是不是我们只须知道如何将训练数据与测试数据输入到模型，然后调用分类或回归的结果就好了呢？
# <p>
# 　　理解其数学原理，自己再实现一遍算法有无必要呢？
# <p>
# 　　在汲取一些前辈的建议以及结合自己的思考后，我觉得自己亲手推一遍数学公式，再用Python实现算法还是有必要的。<br>
# - 因为这样不仅能使我们加深对机器学习本质的理解，还能让我们对模型与数据之间联系更加敏感。<br>
# 例如，如何选择模型、如何选择特征、如何调整参数等等。<br>
# - 其次，目前开源的机器学习算法包提供的都是通用型的算法，并非针对量化投资这一领域来进行优化。<br>
# 所以，当有必要时，我们须根据自己的需求来优化这些通用算法，甚至重写。另一方面，真正有用的算法在量化领域是不会开源的。<br>
# <p>
# 
# 　　当然机器学习这么多的算法也没有必要全都实现一遍，但常见的算法的核心部分还是要亲手实现的。

# ## <b>2、本帖的意义</b>
# ------------------------------------------
# 　　贴主目前只是一个半路入量化坑的在校学生，本身并非研究机器学习这一方向，也不是CS科班出生。<br>
# 但凭借着内心的兴趣还是想扎扎实实地学习机器学习与量化投资的知识，并且希望将来能将两者结合起来做出一些有效的模型或策略。<br>
# 　　本贴是一名普通工科生的自学笔记，<b>笔记详细记录了Logistic Regression的数学原理与推导过程。在核心公式推导过程中，<br>
# 都会引出并详细解释关键的求解技巧和对应的数学知识。因此，只要是对本科高数还有记忆的同学都能完全理解和掌握。</b><br>
# 　　由于贴主才疏学浅，数学、编程、金融样样不通，本帖中错误与不严谨之处在所难免，希望各位朋友批评指正。
# <p>
# 　　最后，本贴剩余部分结构如下：<br>
# 　　3. 介绍Logistic Regression的数学模型，推导并详细解释求解最优回归系数的过程；<br> 
# 　　4. Python实现Logistic Regression的基本版；<br>
# 　　5. 介绍sklearn中的Logistic Regression算法及其关键参数；<br>
# 　　6. 实现一个基于Logistic Regression的简单选股策略。<br>

# ## <b>3、什么是Logistic Regression</b>
# -------------------------------------------------------
# 　　首先，我们知道机器学习是一系列对数据执行分类、回归和聚类等操作的统计算法的统称，这些算法根据历史数据（训练数据）的特征，<br>
# 来对未来数据（测试数据）进行判断。<br>
# 　　本贴我们学习一种最常用也最重要的分类算法—Logistic Regression。<br>
# 　　首先我们引入<b>对数几率函数</b>（logistic function）如下所示：<br>
# $$ y = \frac{1}{1 + e^{-z}} \tag{1} $$
# 　　其函数曲线如下图所示：<br>
#   	![logistic function](https://img3.doubanio.com/view/photo/large/public/p2366159374.jpg) <br>
# 　　从上图我们可以看出，此函数可将定义域为(-∞, +∞)自变量$z$映射到(0,1)区间。若以 $y = 0.5$为阈值，我们可以利用这个函数实现一个二分类器：<br>
# -  当 $z > 0$ ( $y > 0.5$)，将其划分为1类；<br>
# -  当 $z < 0$ ( $y < 0.5$)，将其划分为0类。<br>
# 
# 　　由于于输出$y$为0到1之间的连续值，我们也可以认为函数是对类别的概率估计。
#   <br>
#   <br>
#   <br>
# ### <b>1.1 若有多个输入变量呢？</b>
# 　　那么$z$由下列公式表示：
# $$ z = \omega_0 + \omega_1x_1 + \ldots + \omega_nx_n \tag{2}$$
# 
# 　　其中$\omega_i$表示第$i$个输入变量的权重或是回归系数，第一个输入变量$x_0$默认为1。<br>
#   
# 　　公式(2)的向量形式又如下公式(3)表示：
# $$ z = \boldsymbol{\omega}^\text{T}\boldsymbol{x} \tag{3}$$
# 
# 　　其中，$\boldsymbol{\omega}$与$\boldsymbol{x}$均为$n$行1列的列向量。
# <br>
# <br>
# <br>
# ### <b>1.2 为什么称有如此形式的函数为对数几率函数呢？</b>
# 　　通过简单的推导，Logistic函数可以写成如下形式：
# $$ y = \frac{1}{1+e^{-z}} \Longleftrightarrow \text{ln}\frac{y}{1-y} = \boldsymbol{\omega}^\text{T}\boldsymbol{x} \tag{4}$$
# 　　其中$\text{ln}\frac{y}{1-y}$便称为<b>对数几率</b>（log odds）。“<b>几率</b>”（odds）反映了输入向量$x$被划分为1类的相对可能性。<br>
# 　　因此公式（4）的意义是用线性回归模型去描述类别的<b>对数几率</b>。
# <br>
# <br>
# <br>
# ### <b>1.3 如何求解回归模型中的最优回归系数（权重）呢？</b>
# 　　首先根据公式（4）我们可得：
# $$ \text{ln}\frac{p(y=1|\boldsymbol{x})}{p(y=0|\boldsymbol{x})} = \boldsymbol{\omega}^\text{T}\boldsymbol{x} \tag{5}$$
# 　　进一步我们不难得到在输入为$x$的情况下，输出类别$y$分别为1和0的概率为：
# $$ p(y=1|\boldsymbol{x}) = \frac{e^{ \boldsymbol{\omega}^\text{T}\boldsymbol{x}}}{1 + e^{ \boldsymbol{\omega}^\text{T}\boldsymbol{x}}} = 
# \frac{1}{1 + e^{-\boldsymbol{\omega}^\text{T}\boldsymbol{x}}} \tag{6.1}$$
# $$ p(y=0|\boldsymbol{x}) = \frac{1}{1 + e^{ \boldsymbol{\omega}^\text{T}\boldsymbol{x}}} \tag{6.2} $$
# <br>
# <br>
# <br>
# #### 　　<b>1.3.1 极大似然估计 是求解最优系数的常用方法之一。</b><br>
# 　　极大似然估计的思想为：对于所有的抽样样本，使它们联合概率达到最大的系数便是统计模型最优的系数。<br>
# 　　因此对于第$i$个输入数据$\boldsymbol{x}_i$，它被划分为$y_i$的概率可由公式（7）表示：
# $$ p(y_i|\boldsymbol{x_i};\boldsymbol{\omega}) = p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega})^{y_i}p(y_i=0|\boldsymbol{x_i};\boldsymbol{\omega})^{1-y_i} \tag{7}$$
# 　　公式（7）巧妙之处在于，当$y_i =1$时，$ p(y_i=0|\boldsymbol{x_i};\boldsymbol{\omega})^{1-y_i}=1 $；而当$y_i =0$时, $ p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega})^{y_i}=1 $。<br>
#  <br>
# 　　在极大似然法中，我们假设样本之间都是独立同分布的，因此它们的联合概率就是它们各自概率的乘积。因此关于回归系数$\boldsymbol{\omega}$的极大似然函数可由公式（8）表示：
#   $$  L(\boldsymbol{\omega}) = \prod_{i=1}^{m} p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega})^{y_i}p(y_i=0|\boldsymbol{x_i};\boldsymbol{\omega})^{1-y_i} \tag{8}$$<br>
# 　　为了便于计算，我们对公式（8）两边同时取对数：
# $$  \text{ln}L(\boldsymbol{\omega}) = \text{ln}\left(\prod_{i=1}^{m} p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega})^{y_i}p(y_i=0|\boldsymbol{x_i};\boldsymbol{\omega})^{1-y_i}\right) \tag{9} $$
# 
# 　　根据对数的性质，公式（9）可以逐步写成如下形式：
# \begin{align}
# \text{ln}L(\boldsymbol{\omega}) &= \text{ln}\left(\prod_{i=1}^{m} p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega})^{y_i}p(y_i=0|\boldsymbol{x_i};\boldsymbol{\omega})^{1-y_i}\right)   \nonumber \\
#                                 &= \sum_{i=1}^m y_i\text{ln} p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega}) 
#                                         + (1-y_i)\text{ln} p(y_i=0|\boldsymbol{x_i};\boldsymbol{\omega})       \tag{10}     \\
#                                 &= \sum_{i=1}^m \left( y_i\boldsymbol{\omega}^\text{T}\boldsymbol{x} - \text{ln}(1 + e^{\boldsymbol{\omega}^\text{T}\boldsymbol{x}})\right) \nonumber
# \end{align}
# <br>
# <br>
# <br>
# #### 　　<b>1.3.2 如何求解最优的回归系数 $\boldsymbol{\omega}$ 呢？</b><br>
# 　　有两种经典的数值优化算法：a) 梯度上升法；b) 牛顿法。<br>
# 　　<b>a) 梯度上升法</b>：函数在某一点的梯度总是指向该函数增长最快的方向。因此沿着该函数的梯度方向探寻就能找到该函数的最大值。梯度上升法的迭代公式如下：
# $$ \boldsymbol{\omega}^{t+1} = \boldsymbol{\omega}^{t} + \alpha\nabla_{\boldsymbol{\omega}}\text{ln}L(\boldsymbol{\omega^t}) \tag{11} $$
# 　　公式（11）中 $\alpha$表示每次迭代的步进。$\nabla_{\boldsymbol{\omega}}$表示梯度算子。<br>
# <br>
# 　　根据以上定义，我们求取公式（10）的梯度：<br>
# 　　首先方程两边同时取微分：
# \begin{align}
# d \left( \text{ln}L(\boldsymbol{\omega})\right) 
# &=  \sum_{i=1}^m \left( y_i  d \left(\boldsymbol{\omega}^\text{T}\boldsymbol{x_i} \right) - 
# d \left( \text{ln}(1 + e^{\boldsymbol{\omega}^\text{T}\boldsymbol{x_i}})\right) \right)           \nonumber         \\
# &= \sum_{i=1}^m \left( y_i  d \left( \boldsymbol{\omega}^\text{T}\boldsymbol{x_i} \right)^\text{T} - 
# \frac{e^{ \boldsymbol{\omega}^\text{T}\boldsymbol{x_i}}}{1 + e^{ \boldsymbol{\omega}^\text{T}\boldsymbol{x_i}}}
# d \left( \boldsymbol{\omega}^\text{T}\boldsymbol{x_i} \right)^\text{T} \right)                    \nonumber           \\
# &= \sum_{i=1}^m y_i\boldsymbol{x_i}^\text{T}d(\boldsymbol{\omega}) -
# p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega})\boldsymbol{x_i}^\text{T}d(\boldsymbol{\omega})       \nonumber           \\
# &= \sum_{i=1}^m \left( \boldsymbol{x_i}^\text{T}(y_i - p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega}))   \right) d(\boldsymbol{\omega})   \tag{12}
# \end{align}
# 
# 　　因此公式（10）的梯度最终可写成公式（13）的形式：
# $$  
# \nabla_\boldsymbol{\omega}\text{ln}L(\boldsymbol{\omega}) = 
# \frac{\partial\text{ln}L(\boldsymbol{\omega})}{\partial \boldsymbol{\omega}} = 
# \sum_{i=1}^m  \boldsymbol{x_i}(y_i - p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega}))     \tag{13}
# $$
# <br>
# 　　从公式（12）到（13）有两处向量微分运算的小技巧（刚开始我也不会...）：<br>
# 　　第一处：公式（12）的第二行，由于$\boldsymbol{\omega}^\text{T}\boldsymbol{x_i}$为标量，我们取转置后，其性质不变。<br>
# 　　第二处：在向量的微分运算中，梯度与导数矩阵是相互装置的关系（这样说可能不是特别严谨...）。<br>
# 　　　　　　所以公式（12）的$\boldsymbol{x_i}^\text{T}$在公式（13）中变为$\boldsymbol{x_i}$。<br>
# <br>
# <br>
# 　　<b>b) 牛顿法</b>：其原理是利用泰勒公式不断迭代，从而逐次逼近零点或极值点。<br>
# 　　<b>一阶展开求零点</b>：<br>
# 　　我们知道泰勒公式在$x_0$处的一阶展开如下所示：<br>
# $$ f(x) \thickapprox f(x_0) + (x - x_0)f^{'}(x_0)  \tag{14} $$
# 　　我们用一阶展开式求解函数的零点：
# $$ f(x_0) + (x - x_0)f^{'}(x_0) = 0  \tag{15} $$
# 　　稍作整理，公式（15）可写成如下形式：
# $$ x = x_0 - \frac{f(x_0)}{f^{'}(x_0)} \tag{16}$$
# 　　由于一阶展开式只是与原函数近似相等，因此公式（16）中的$x$并非函数的零点，而只是比$x_0$更接近零点。因此通过对公式（16）迭代可以不断逼近函数零点。<br>
# <br>
# <br>
# 　　<b>如果是求函数的极值点呢？</b><br>
# 　　那我们就要用到泰勒二阶展开式：
# $$ f(x) \thickapprox f(x_0) + (x - x_0)f^{'}(x_0) + \frac{1}{2}(x - x_0)^2f^{"}(x_0) = g(x)\tag{17} $$
# 　　我们知道函数在其一阶导数等于0时取到极值，因此我们求上式(17) $ g^{'}(x) = 0$的解：
# $$ g^{'}(x) = f^{'}(x_0) + (x - x_0)f^{"}(x_0) = 0 \tag{18} $$
# 　　通过像公式（15）一样整理上式（18），可得：
#   $$ x = x_0 - \frac{f^{'}(x_0)}{f^{"}(x_0)} \tag{19}$$ <br>
# <br>  
# <br>
# 　　<b>如果输入是多维的变量呢？</b><br>
# 　　高维情况的牛顿法迭代公式与公式（19）相似，其如下所示：
# $$ \boldsymbol{x}^{t+1} = \boldsymbol{x}^{t} - \left[ Hf(\boldsymbol{x}^{t}) \right]^{-1}\nabla f(\boldsymbol{x}^{t}) \tag{20}$$
# 　　其中$\nabla f(\boldsymbol{x}^{t})$表示的是函数$f(\boldsymbol{x}^{t})$的梯度，$\left[ Hf(\boldsymbol{x}^{t}) \right]^{-1}$表示的是函数的	Hessian矩阵的逆矩阵。<br>
# 　　当输入变量为一维时，Hessian矩阵其实可以理解为变量的二阶导数。其Hessian矩阵的定义如下：
# $$
# H(f(\boldsymbol{x})) = 
# \frac{\partial^{2}f(\boldsymbol{x})}{\partial{\boldsymbol{x}}\partial{\boldsymbol{x}}^{\text{T}}} = 
# \left[
# \begin{matrix}
# \frac{\partial^{2}f}{\partial{x_1^2}}  & \frac{\partial^{2}f}{\partial{x_1}\partial{x_2}}  & \cdots & \frac{\partial^{2}f}{\partial{x_1}\partial{x_n}}    \\
# \frac{\partial^{2}f}{\partial{x_2}\partial{x_1}}  & \frac{\partial^{2}f}{\partial{x_2^2}}  & \cdots & \frac{\partial^{2}f}{\partial{x_2}\partial{x_n}}    \\
# \vdots & \vdots & \ddots & \vdots     \\
# \frac{\partial^{2}f}{\partial{x_n}\partial{x_1}}  & \frac{\partial^{2}f}{\partial{x_n}\partial{x_2}} & \cdots & \frac{\partial^{2}f}{\partial{x_n^2}}     \\
# \end{matrix}
# \right]
# \tag{21}$$
# <br>
# <br>
# <br>
# 　　<b>梯度上升法与牛顿法的比较：</b><br>
# 　　梯度法是一阶收敛的，而牛顿法是二阶收敛的。因此牛顿法的迭代次数要少于梯度法，然而Hessian的逆矩阵的计算会增加算法的复杂度。<br>
# 　　这个问题又可以通过Quasi-Newton方法解决，本帖对此方法也不做深入讨论。<br>
# <br>

# ## <b>4. Python实现Logistic Regression的基本版</b>
# -------------------------------------------------------
# 　　我们选择牛顿法求解Logistic Regression中的最优回归系数。<br>
# 　　根据牛顿法迭代公式（20）和Hessian矩阵的定义（21），本贴Logistic Regression中的回归系数迭代更新公式为：<br>
# $$
# \boldsymbol{\omega}^{t+1} = \boldsymbol{\omega}^{t} - 
# \left[ 
# \frac{\partial^2\text{ln}L(\boldsymbol\omega)} { \partial\boldsymbol\omega\partial\boldsymbol\omega^{\text{T}}}
# \right]^{-1}
# \frac{\partial\text{ln}L(\boldsymbol\omega)}{\partial\boldsymbol\omega}
# \tag{22}
# $$
# 　　其中$\frac{\partial\text{ln}L(\boldsymbol\omega)}{\partial\boldsymbol\omega}$在公式（13）中已求得，而对$\frac{\partial^2\text{ln}L(\boldsymbol\omega)} { \partial\boldsymbol\omega\partial\boldsymbol\omega^{\text{T}}}$的求解可按照上文所提到的两条向量微分运算方法来计算，因此：
# $$
# \frac{\partial^2\text{ln}L(\boldsymbol\omega)} { \partial\boldsymbol\omega\partial\boldsymbol\omega^{\text{T}}} = 
# \sum_{i=1}^m \boldsymbol x_i \boldsymbol x_i^{\text{T}} p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega}) \left( p(y_i=1|\boldsymbol{x_i};\boldsymbol{\omega}) - 1 \right)
# \tag{23}
# $$
# 　　好了！到此迭代公式已经完全掌握了，现在该实现算法了！<br>
# <br>
# <br>
# <br>
# ### <b>4.1 准备数据：</b>
# 　　从网上整理了一组数据: 前两列元素是feature；最后一列是label，其值是0或1。<br>
# ```Python
# data = [[6.9, 3.1, 1.0], [5.0, 3.5, 0.0], [5.0, 2.3, 1.0], [4.6, 3.4, 0.0], [5.5, 2.4, 1.0], [5.5, 3.5, 0.0], [5.2, 4.1, 0.0], [4.7, 3.2, 0.0], [5.4, 3.9, 0.0], [5.6, 3.0, 1.0], [6.1, 2.8, 1.0], [5.8, 2.7, 1.0], [6.2, 2.2, 1.0], [5.0, 3.4, 0.0], [5.7, 3.8, 0.0], [6.3, 2.5, 1.0], [6.0, 3.4, 1.0], [4.9, 3.1, 0.0], [4.9, 3.0, 0.0], [5.6, 2.7, 1.0], [5.7, 2.8, 1.0], [5.0, 3.2, 0.0], [5.9, 3.0, 1.0], [6.1, 2.8, 1.0], [5.1, 3.5, 0.0], [4.6, 3.2, 0.0], [5.1, 2.5, 1.0], [5.1, 3.7, 0.0], [5.4, 3.7, 0.0], [5.3, 3.7, 0.0], [4.8, 3.0, 0.0], [4.3, 3.0, 0.0], [4.9, 2.4, 1.0], [5.5, 2.5, 1.0], [5.8, 4.0, 0.0], [5.1, 3.8, 0.0], [5.7, 2.8, 1.0], [5.1, 3.8, 0.0], [5.0, 3.0, 0.0], [4.4, 3.0, 0.0], [5.4, 3.4, 0.0], [5.7, 2.6, 1.0], [6.4, 3.2, 1.0], [5.2, 3.5, 0.0], [4.8, 3.4, 0.0], [6.0, 2.7, 1.0], [5.6, 3.0, 1.0], [5.1, 3.3, 0.0], [5.4, 3.9, 0.0], [6.8, 2.8, 1.0], [5.9, 3.2, 1.0], [4.5, 2.3, 0.0], [6.1, 3.0, 1.0], [4.4, 2.9, 0.0], [5.4, 3.4, 0.0], [5.5, 2.4, 1.0], [5.6, 2.9, 1.0], [5.5, 2.6, 1.0], [5.7, 4.4, 0.0], [5.0, 3.6, 0.0], [5.0, 3.4, 0.0], [4.6, 3.6, 0.0], [5.7, 2.9, 1.0], [5.0, 2.0, 1.0], [6.3, 2.3, 1.0], [7.0, 3.2, 1.0], [4.6, 3.1, 0.0], [5.0, 3.3, 0.0], [5.1, 3.8, 0.0], [4.9, 3.1, 0.0], [4.7, 3.2, 0.0], [5.5, 2.3, 1.0], [6.7, 3.1, 1.0], [6.5, 2.8, 1.0], [6.7, 3.0, 1.0], [4.4, 3.2, 0.0], [6.6, 2.9, 1.0], [5.0, 3.5, 0.0], [6.4, 2.9, 1.0], [5.1, 3.5, 0.0], [6.3, 3.3, 1.0], [5.4, 3.0, 1.0], [6.2, 2.9, 1.0], [5.1, 3.4, 0.0], [5.2, 3.4, 0.0], [4.8, 3.0, 0.0], [6.0, 2.9, 1.0], [5.5, 4.2, 0.0], [4.8, 3.4, 0.0], [5.8, 2.6, 1.0], [5.8, 2.7, 1.0], [6.7, 3.1, 1.0], [4.8, 3.1, 0.0], [5.7, 3.0, 1.0], [6.1, 2.9, 1.0], [6.0, 2.2, 1.0], [6.6, 3.0, 1.0], [4.9, 3.1, 0.0], [5.6, 2.5, 1.0], [5.2, 2.7, 1.0]]
# ```
# <br>
# ### <b>4.2 数据预处理，包括：</b>
# 　　1.增加常数项$x_0$: 1.0；<br>
# 　　2.feature 与 label的分离；<br>
# 　　3.返回数组类型的feature 和 label。
# ```Python
# def preprocess(data):
#     # labels 用来标识每个数据的类别:0或1
#     labels = [] 
#     # train_X 输入数据列表: [[1.0, x_i1, x_i2],..., [1.0, x_n1, x_n2]]
#     train_X = [] 
#     for x in data:
#         # 二维平面的直线方程为： ax + by + c = 0，因此我们需要再原来的输入数据中增加一常数项 1.0
#         train_X.append([1.0, x[0], x[1]])
#         labels.append(x[2])
#     return np.array(train_X), np.array(labels)
# ```
# <br>
# ### <b>4.3 定义Logistic函数：</b>
# ```Python
# def logit(z):
#     '''
#         logistic function 即：
#         牛顿迭代公式中的p(yi=1 | xi; w)
#     '''
#     return 1.0 / (1.0 + np.exp(-z))
# ```
# <br>
# ### <b>4.4 梯度上升法算法：</b>
# ```Python
# def grad(train_X, labels, max_iter = 500):
# 	# 下面两行代码用于矩阵化 训练数据
#     X = np.mat(train_X) # 100 行 3列
#     Y = np.mat(labels).transpose() # 100 行 1列
#     rows, columns = np.shape(X)
#     # 初始回归系数 weights 为全1的向量, weights的行数与输入的训练数据维度相同
#     weights = np.ones((columns, 1))
#     # 步进alpha: 0.001
#     alpha = 0.001
#     # t 表示当前迭代的次数
#     t = 0
#     while t < max_iter:
#         P1 = logit(np.dot(X, weights))
#         weights += alpha * X.T * (Y - P1)
#         t += 1
#     
#     return weights
# ```
# 　　我们还是给出了梯度上升法的算法实现。<br>
# 　　上述算法的倒数第三行，我们用矩阵形式表示梯度上升法的迭代公式，其具体数学形式如下：
# $$
# \left[
# \begin{matrix}
# \omega_0  \\
# \omega_1  \\
# \omega_2 
# \end{matrix}
# \right]^{(t+1)} 
# = 
# \left[
# \begin{matrix}
# \omega_0  \\
# \omega_1  \\
# \omega_2 
# \end{matrix}
# \right]^{(t)} 
# + 
# \alpha
# \left[
# \begin{matrix}
# 1.0 & x_{11} & x_{12} \\
# 1.0 & x_{21} & x_{22} \\
# \vdots & \vdots & \vdots \\
# 1.0 & x_{n1} & x_{n2} \\
# \end{matrix}
# \right]^{\text{T}} 
# 
# \left[ 
# \left(
# \begin{matrix}
# y_1 \\
# y_2 \\
# \vdots \\
# y_n   \\
# \end{matrix}
# \right) 
# -
# \left( 
# \begin{matrix}
# p_1(\boldsymbol{x}_1; \boldsymbol\omega) \\
# p_1(\boldsymbol{x}_2; \boldsymbol\omega) \\
# \vdots \\
# p_1(\boldsymbol{x}_n; \boldsymbol\omega)
# \end{matrix}
# \right)
# \right]
# $$
# <br>
# ### <b>4.5 牛顿法算法：</b>
# ```Python
# def newton(train_X, labels, max_iter = 10):
#     X = np.mat(train_X)
#     Y = np.mat(labels).transpose()
#     rows, columns = np.shape(X)
#     # 初始回归系数 weights 为全0的向量, weights的行数与输入的训练数据维度相同。
#     weights = np.zeros((columns, 1))
#     # t 表示当前迭代的次数
#     t = 0
#     while t < max_iter:
#         P1 = logit(np.dot(X, weights))
#         gradient =  np.dot(X.T, (Y - P1))
#         P_mat = np.array(P1) * np.array(P1 - 1) * np.eye(rows)
#         hessian = X.T * P_mat * X
#         # np.linalg.inv() 是求的Hessian矩阵的逆矩阵
#         weights -= np.linalg.inv(hessian) * gradient
#         t += 1
#         
#     return weights
# ```
# 　　与上小节4.4一致，我们都用矩阵形式表示迭代公式。其中梯度的具体数学形式在上小节已给出。<br>
# 　　在本节的牛顿算法中，Hessian矩阵具体数学形式如下：
# $$
# Hessian = 
# \left[
# \begin{matrix}
# 1.0 & x_{11} & x_{12} \\
# 1.0 & x_{21} & x_{22} \\
# \vdots & \vdots & \vdots \\
# 1.0 & x_{n1} & x_{n2} \\
# \end{matrix}
# \right]^{\text{T}} 
# 
# \left[
# \begin{matrix}
# p_1(\boldsymbol{x}_1; \boldsymbol\omega) (p_1(\boldsymbol{x}_1; \boldsymbol\omega) - 1)  & \\
# & p_1(\boldsymbol{x}_2; \boldsymbol\omega) (p_1(\boldsymbol{x}_2; \boldsymbol\omega) - 1) &\\
# & \ddots \\
# & &  p_1(\boldsymbol{x}_n; \boldsymbol\omega) (p_1(\boldsymbol{x}_n; \boldsymbol\omega) - 1)
# \end{matrix}
# \right]
# 
# \left[
# \begin{matrix}
# 1.0 & x_{11} & x_{12} \\
# 1.0 & x_{21} & x_{22} \\
# \vdots & \vdots & \vdots \\
# 1.0 & x_{n1} & x_{n2} \\
# \end{matrix}
# \right]
# $$
# <br>
# ### <b>4.6 选择最优算法，根据求得的最优系数画出回归直线：</b>
# ```Python
# # 牛顿法
# weights = newton(preprocess(data)[0], preprocess(data)[1])
# # 梯度法
# #weights = grad(preprocess(data)[0], preprocess(data)[1])
# # 根据数据的类别，得到两个数据集
# x0=list()
# y0=list()
# x1=list()
# y1=list()
# for i in data:
#     if i[-1] == 0:
#         x0.append(i[0])
#         y0.append(i[1])
#     else:
#         x1.append(i[0])
#         y1.append(i[1])
# 
# # c=sns.color_palette()[1] ： seaborn 包中的配色方案
# plt.scatter(x0, y0, marker='^', c=sns.color_palette()[1])
# plt.scatter(x1, y1, marker='o', c=sns.color_palette()[2])
# 
# x = np.arange(4.0, 7.0, 0.1)
# # 画出回归直线： w0 + w1*x1 + w2*x2 = 0, 我们定义x1为横轴坐标，x2为纵轴坐标
# y = (-weights[0] - weights[1]*x)/weights[2]
# plt.plot(x, y, color=sns.color_palette()[0])
# plt.xlabel('x1')
# plt.ylabel('x2')
# plt.show()
# ```
# 
# ### <b>4.7 下面一个cell运行整段代码：</b>

# In[ ]:

import numpy as np #数值计算包
import matplotlib.pyplot as plt #绘图包
import seaborn as sns # 绘图美化包 用于配色

data = [[6.9, 3.1, 1.0], [5.0, 3.5, 0.0], [5.0, 2.3, 1.0], [4.6, 3.4, 0.0], [5.5, 2.4, 1.0], [5.5, 3.5, 0.0], [5.2, 4.1, 0.0], [4.7, 3.2, 0.0], [5.4, 3.9, 0.0], [5.6, 3.0, 1.0], [6.1, 2.8, 1.0], [5.8, 2.7, 1.0], [6.2, 2.2, 1.0], [5.0, 3.4, 0.0], [5.7, 3.8, 0.0], [6.3, 2.5, 1.0], [6.0, 3.4, 1.0], [4.9, 3.1, 0.0], [4.9, 3.0, 0.0], [5.6, 2.7, 1.0], [5.7, 2.8, 1.0], [5.0, 3.2, 0.0], [5.9, 3.0, 1.0], [6.1, 2.8, 1.0], [5.1, 3.5, 0.0], [4.6, 3.2, 0.0], [5.1, 2.5, 1.0], [5.1, 3.7, 0.0], [5.4, 3.7, 0.0], [5.3, 3.7, 0.0], [4.8, 3.0, 0.0], [4.3, 3.0, 0.0], [4.9, 2.4, 1.0], [5.5, 2.5, 1.0], [5.8, 4.0, 0.0], [5.1, 3.8, 0.0], [5.7, 2.8, 1.0], [5.1, 3.8, 0.0], [5.0, 3.0, 0.0], [4.4, 3.0, 0.0], [5.4, 3.4, 0.0], [5.7, 2.6, 1.0], [6.4, 3.2, 1.0], [5.2, 3.5, 0.0], [4.8, 3.4, 0.0], [6.0, 2.7, 1.0], [5.6, 3.0, 1.0], [5.1, 3.3, 0.0], [5.4, 3.9, 0.0], [6.8, 2.8, 1.0], [5.9, 3.2, 1.0], [4.5, 2.3, 0.0], [6.1, 3.0, 1.0], [4.4, 2.9, 0.0], [5.4, 3.4, 0.0], [5.5, 2.4, 1.0], [5.6, 2.9, 1.0], [5.5, 2.6, 1.0], [5.7, 4.4, 0.0], [5.0, 3.6, 0.0], [5.0, 3.4, 0.0], [4.6, 3.6, 0.0], [5.7, 2.9, 1.0], [5.0, 2.0, 1.0], [6.3, 2.3, 1.0], [7.0, 3.2, 1.0], [4.6, 3.1, 0.0], [5.0, 3.3, 0.0], [5.1, 3.8, 0.0], [4.9, 3.1, 0.0], [4.7, 3.2, 0.0], [5.5, 2.3, 1.0], [6.7, 3.1, 1.0], [6.5, 2.8, 1.0], [6.7, 3.0, 1.0], [4.4, 3.2, 0.0], [6.6, 2.9, 1.0], [5.0, 3.5, 0.0], [6.4, 2.9, 1.0], [5.1, 3.5, 0.0], [6.3, 3.3, 1.0], [5.4, 3.0, 1.0], [6.2, 2.9, 1.0], [5.1, 3.4, 0.0], [5.2, 3.4, 0.0], [4.8, 3.0, 0.0], [6.0, 2.9, 1.0], [5.5, 4.2, 0.0], [4.8, 3.4, 0.0], [5.8, 2.6, 1.0], [5.8, 2.7, 1.0], [6.7, 3.1, 1.0], [4.8, 3.1, 0.0], [5.7, 3.0, 1.0], [6.1, 2.9, 1.0], [6.0, 2.2, 1.0], [6.6, 3.0, 1.0], [4.9, 3.1, 0.0], [5.6, 2.5, 1.0], [5.2, 2.7, 1.0]]

def preprocess(data):
    # labels 用来标识每个数据的类别:0或1
    labels = [] 
    # train_X 输入数据列表: [[1.0, x_i1, x_i2],..., [1.0, x_n1, x_n2]]
    train_X = [] 
    for x in data:
        # 二维平面的直线方程为： ax + by + c = 0，因此我们需要再原来的输入数据中增加一常数项 1.0
        train_X.append([1.0, x[0], x[1]])
        labels.append(x[2])
    return np.array(train_X), np.array(labels)


def logit(z):
    '''
        logistic function 即：
        牛顿迭代公式中的p(yi=1 | xi; w)
    '''
    return 1.0 / (1.0 + np.exp(-z))
    
def grad(train_X, labels, max_iter = 500):
    X = np.mat(train_X) # 100 行 3列
    Y = np.mat(labels).transpose() # 100 行 1列
    rows, columns = np.shape(X)
    # 初始回归系数 weights 为全1的向量, weights的行数与输入的训练数据维度相同
    weights = np.ones((columns, 1))
    # 步进alpha: 0.001
    alpha = 0.001
    # t 表示当前迭代的次数
    t = 0
    while t < max_iter:
        P1 = logit(np.dot(X, weights))
        weights += alpha * X.T * (Y - P1)
        t += 1
    
    return weights

def newton(train_X, labels, max_iter = 10):
    X = np.mat(train_X)
    Y = np.mat(labels).transpose()
    rows, columns = np.shape(X)
    # 初始回归系数 weights 为全0的向量, weights的行数与输入的训练数据维度相同。
    weights = np.zeros((columns, 1))
    # t 表示当前迭代的次数
    t = 0
    while t < max_iter:
        P1 = logit(np.dot(X, weights))
        gradient =  np.dot(X.T, (Y - P1))
        P_mat = np.array(P1) * np.array(P1 - 1) * np.eye(rows)
        hessian = X.T * P_mat * X
        weights -= np.linalg.inv(hessian) * gradient
        t += 1
    
    return weights
    
# 牛顿法
weights = newton(preprocess(data)[0], preprocess(data)[1])
# 梯度法
# weights = grad(preprocess(data)[0], preprocess(data)[1])
# 根据数据的类别，得到两个数据集
x0=list()
y0=list()
x1=list()
y1=list()
for i in data:
    if i[-1] == 0:
        x0.append(i[0])
        y0.append(i[1])
    else:
        x1.append(i[0])
        y1.append(i[1])

# c=sns.color_palette()[1] ： seaborn 包中的配色方案
plt.scatter(x0, y0, marker='^', c=sns.color_palette()[1])
plt.scatter(x1, y1, marker='o', c=sns.color_palette()[2])

x = np.arange(4.0, 7.0, 0.1)
# 画出回归直线： w0 + w1*x1 + w2*x2 = 0, 我们定义x1为横轴坐标，x2为纵轴坐标
y = (-weights[0] - weights[1]*x)/weights[2]
plt.plot(x, y, color=sns.color_palette()[0])
plt.xlabel('x1')
plt.ylabel('x2')
plt.show() 


# ## <b>5. sklearn中的Logistic Regression算法：</b>
# ---------------------------------------
# 　　sklearn是Python实现的开源机器学习算法包，优矿正好也支持它。下面简单介绍下sklearn中Logistic Regression的使用。<br>
# <br>
# ### <b>5.1 实例化Logistic Regression：</b>
# ```Python 
# from sklearn.linear_model import LogisticRegression
# # 参数C为正则化强度的倒数，用于控制回归系数的复杂度。
# lr = LogisticRegression(C = 1.0)
# ```
# 　　理解什么是<b>正则化</b>，我们先理解正则化的作用：<b> 正则化可防止模型过于复杂。</b><br>
# 　　我们知道机器学习算法本质是一个最优化问题，算法根据训练数据所求得的模型系数可以使得数据分类、拟合等准确率更高。<br>
# 　　在保证模型准确性的同时，我们同样希望模型不过于复杂，因此我们就要在两者之间权衡。<br>
# 　　模型的准确性一般用代价函数来衡量，而模型的复杂度一般用系数的平方表示（系数向量的内积）：
# $$ L^{'}(\boldsymbol\omega) = L(\boldsymbol\omega) + C \boldsymbol\omega^\text{T}\boldsymbol\omega \tag{24}$$ 
# 　　在本问题中，我们将极大似然函数的负值视为代价函数，代价函数越小，说明极大似然函数越大，说明模型分类的准确性越高；<br>
# 　　$C\boldsymbol\omega^\text{T}\boldsymbol\omega$ 用以衡量模型复杂度，$C$用以控制原代价函数与模型复杂度之间的折中。<br>
# 　　<b>很显然，$C$越大，算法对模型复杂度越敏感。</b>
# <br>
# <br>
# ### <b>5.2 训练数据：</b>
# ```Python
# lr.fit(train_X, labels)
# ```
# 　　其中，train_X 表示的训练数据的特征，其格式$n$行$m$列的矩阵，$n$为训练数据的个数，$m$为数据的特征个数（维度）；<br>
# 　　labels 表示的训练数据的类别，格式为$n$行1列的列向量。
# <br>
# <br>
# ### <b>5.3 输出模型系数，输出测试数据的类别与概率：</b>
# ```Python
# # 回归系数
# print '输出模型回归系数：', lr.coef_
# print '\n'
# 
# # 预测类别
# print '输入数据[1.0, 5.0, 2.0]的类别：', lr.predict([1.0, 5.0, 2.0])
# print '输入数据[1.0, 5.0, 4.0]的类别：', lr.predict([1.0, 5.0, 4.0])
# print '\n'
# 
# # 预测概率
# print  '输入数据[1.0, 5.0, 2.0]为0类1类的概率：', lr.predict_proba([1.0, 5.0, 2.0])
# print  '输入数据[1.0, 5.0, 4.0]为0类1类的概率：', lr.predict_proba([1.0, 5.0, 4.0])
# ```
# 　　如果用上一节准备好数据，我们可分别得到如下结果：
# ##### <b>回归系数</b>
# ```Python
# 输出模型回归系数： [[-0.55312625  2.2765652  -3.63240141]]
# 
# ```
# 
# ##### <b>预测类别</b>
# ```Python
# 输入数据[1.0, 5.0, 2.0]的类别： [ 1.]
# 输入数据[1.0, 5.0, 4.0]的类别： [ 0.]
# 
# ```
# 
# ##### <b>预测概率</b>
# ```Python
# 输入数据[1.0, 5.0, 2.0]为0类1类的概率： [[ 0.04689694  0.95310306]]
# 输入数据[1.0, 5.0, 4.0]为0类1类的概率： [[ 0.98597835  0.01402165]]
# 
# ```
# 

# ## <b>6. 基于Logistic Regression的简单选股策略：</b>
# ---------------------------------------
# 　　该策略原理简单，即用T-1天股票的因子值，预测10个交易日后股票的涨跌。<br>
# 　　训练数据为T-1 到 T-60 之间的滚动时间窗口。<br>
# 　　策略示意图如下所示：<br>
# $$  \underbrace{T-60,  \dots, \overbrace{\overbrace{T-11 ,}^{\mbox{取特征因子}} \dots\dots,\overbrace{T-1,}^{\mbox{判断类别：1涨0跌}}}^{\mbox{一组训练数据}}}_{\mbox{所有训练数据}}     T   $$
# <br>
# 　　在回测时，为了避免幸存者偏差，选股都是使用set_universe('HS300', begin_data)，没有直接使用account.universe。<br>
# <br>
# 　　下面的cell给出策略代码及回测结果：
